function [correl_dX_Y,covar_dX_Y,std_dX,std_Y] = CalcTimeSeriesCorr3(X,Y,prob_X,Phi_X,n,Phi_X_n)
% computes corr(X(t+m)-X(t),Y(t+m))
    
    if nargin<=5
        Phi_X_n = mpower2(Phi_X,n);
    end
    
    mean_X = sum(prob_X(:).*X(:));
    var_X = sum(prob_X(:).*(X(:).^2)) - mean_X^2;
    mean_Y = sum(prob_X(:).*Y(:));
    std_Y = sqrt(sum(prob_X(:).*(Y(:).^2)) - mean_Y^2);
    
    condE_X_n = Phi_X_n*X(:);
    var_dX = 2*var_X - 2*(sum(prob_X(:).*X(:).*condE_X_n) - mean_X^2);
    std_dX = sqrt(var_dX);
    
    condE_Y_n = Phi_X_n*Y(:);
    
    covar_dX_Y = sum(prob_X(:).*(X(:).*Y(:) - X(:).*condE_Y_n(:)));
    
    correl_dX_Y = covar_dX_Y/(std_dX*std_Y);


end